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Abstract 



■ We introduce a numerical method for solving Grad's moment equations or regular- 

ized moment equations for arbitrary order of moments. In our algorithm, we do not 
explicitly need the moment equations. Instead, we directly start from the Boltzmann 
equation and perform Grad's moment method [12] and the regularization technique 
[27] numerically. We dchne a conservative projection operator and propose a fast 
implementation which makes it convenient to add up two distributions and provides 
more efficient flux calculations compared with the classic method using explicit ex- 
pressions of flux functions. For the collision term, the BGK model is adopted so that 
the production step can be done trivially based on the Hermite expansion. Extensive 
numerical examples for one- and two-dimensional problems are presented. Conver- 
gence in moments can be validated by the numerical results for different number of 
—i ' moments. 
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1 Introduction 



In recent years, the simulation of rarefied fluids or microflows, which contain significant 
non-equilibrium characteristics, became one of the major directions of fluid dynamics. 
The Boltzmann equation, which is considered as the basis of modern kinetic theory, is the 
starting point of such simulations. Because of the high dimension of variables and the 
complicated form of its collision operators, people tend to use its discrete or simplified 
form instead of the Boltzmann equation itself in numerical simulation. Lots of work has 
been done to simplify the collision operator, such as the BGK model [4], the Shakhov 
model [23], the ES-BGK model [17], the Liu's model [18], the Maxwell molecules model 
[11], and so on. Another way of simplification is to discretize the Boltzmann equation 
by some expansion. In this field, the Chapman-Enskog expansion [8, 10] and the Grad's 
expansion [12, 14] achieved great success in the early exploration of kinetic theory. 

However, both methods of expansion suffer some problems which greatly restrict their 
application. The Chapman-Enskog expansion shows an unstable behavior in the case 
of high order expansions, such as Burnett and super-Burnett equations [5]; the Grad's 
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moment equations lead to unphysical subshocks when the Mach number is large (see 
e.g. [26]). In order to extend their use, several corrections are applied to these models, 
which include the R13 model, proposed by H. Struchtrup and M. Torrilhon in [27]. The 
regularization of Grad's moment equations is done by combining Grad's technique with 
first-order Chapman-Enskog expansion. [38] summarizes the recent work on R13 equa- 
tions, in which it is mentioned that such a regularization technique can be extended to 
any system obtained by Grad's moment method. In this paper, we introduce a method 
that numerically solves the Grad's equations for arbitrary order of moments together with 
their regularization. The numerical method generating the Grad's moment equations for 
arbitrary number of moments has been proposed in [34] and implemented in [2]. And the 
results for shock tube are reported in [33, 3]. However, we solve these equations "only 
numerically," which means we do not need the explicit forms of those equations in our 
algorithm. The regularized moment equations are also considered, and to our knowledge, 
no results about the generation of regularized equations for arbitrary number of moments 
have been reported. Numerical methods for R13 equations are discussed in [32, 15, 37], 
and the computational framework for the R20 equations can be found in [21]. 

Our starting point is the Boltzmann equation, and we adopt the BGK model and a uni- 
form rectangular mesh for discretizing the spatial variable for simplicity. The distribution 
function is expanded in an Hermite series using the method in [12]. The series is truncated 
at a certain place and the coefficients in the expansion are stored for each cell. Following 
the standard procedure, we adopt the classic time splitting method in our scheme. The 
convection term of the Boltzmann equation is discretized by the finite volume method 
with the HLL numerical fluxes. This requires algebraic operations on the distributions 
between the neighbouring cells. Grad's method expands the distribution function by the 
Hermite functions with center at the local mean velocity and a scaling factor associated 
with the local temperature, while these parameters are different in different places. This 
makes it nontrivial to add up two distributions in different cells. Thus it is extremely com- 
plicated to calculate the numerical fluxes directly. As a key technique in this paper, we 
propose a fast algorithm which projects a distribution function expanded in the discrete 
space with one mean velocity and macroscopic temperature to that with another mean 
velocity and macroscopic temperature. This projection is conservative with respect to all 
the moments that are not truncated. Moreover, we prove the projection is invertible so 
that no information will be lost. With the help of this fast projection, distributions in 
any two neighbouring cells can be transformed into the same space efficiently. In order to 
estimate the signal velocities in the HLL numerical fluxes, we prove that the eigenvalues 
of the Jacobian matrix of the linearized flux function are actually the roots of Hermite 
functions plus the mean velocity. Thus we can take the eigenvalues with maximal absolute 
values as the approximated signal velocities. Now, standard HLL numerical fluxes can be 
calculated conveniently. After accumulating the contribution of the convection term, the 
expansion is not of Grad's type any more, since Grad's expansion requires all first-order 
coefficients and the trace of second-order coefficients to be zero. Again, the fast projection 
is applied to correct the center and scaling factor of the expansion so that the properties of 
Grad's expansion can be recovered. These properties make it trivial to perform production 
of the BGK model by a direct scaling of the moments with orders not lower than two. 

In [27], H. Struchtrup and M. Torrilhon used the Chapman-Enskog expansion to de- 
duce the regularized system of Grad's 13-moment equations — R13 equations. The basic 
idea of this regularization can be viewed as a strategy to "guess" the truncated moments. 
Using this idea, we apply the Chapman-Enskog expansion of the Boltzmann equation 
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around the Grad non-equilibrium manifolds [26] of arbitrary order. The closure of the 
system is then achieved by the standard asymptotic techniques therein. This method 
can be perfectly integrated into our numerical scheme without deducing the macroscopic 
equations by intricate algebraic calculations. For arbitrary order of moments, the reg- 
ularization in our algorithm introduces only first order derivative terms, which can be 
numerically approximated by gradient reconstruction. 

In our method, the computational cost of the fast conservative projection is linear in 
terms of the number of moments. Thus it is essentially faster to calculate the numerical 
fluxes than the classic method using the flux functions in the macroscopic equations. 
Actually, the macroscopic equations have never been deduced, but they have been solved 
implicitly by our method. Then the framework of our method can appear to be uniform 
for moment equations of any order. This makes it very convenient to implement our 
algorithm. We need not deduce and code for the complicated flux functions at all in the 
case of high order. Moreover, we need only to develop one copy of the code for all different 
orders. 

We carry out numerical experiments in both one- and two-dimensional cases. Different 
Knudsen numbers and different orders of moments are examined to demonstrate the use- 
fulness of large moment systems, and the convergence in moments is validated numerically. 
Regularized moment systems ranging from 20 moments up to 455 moments are simulated 
in our one-dimensional examples. Two-dimensional examples for up to 84 moments are 
presented, and there is even an example demonstrating the capacity of our method to sim- 
ulate a three-dimensional non-equilibrium process. To the best of our knowledge, it is the 
first time that the method for arbitrary order regularized moment equations is numerically 
implemented, and the moment method for large systems are applied to two-dimensional 
problems. 

The layout of this paper is as follows: in Section 2, an overview of Boltzmann equation 
and the BGK model is given as the basis of our algorithm. In section 3, the details of the 
algorithm which generate numerical solution for arbitrary order Grad's moment equations 
are introduced. In section 4, regularization of moment equations is considered and then 
the whole algorithm is outlined. We present in Section 5 four numerical examples including 
one- and two-dimensional tests to make a comparison between results for different moment 
equations, different Knudsen numbers and different meshes. At last, some concluding 
remarks will be given in Section 6. 

2 The Boltzmann equation and BGK collision model 

In the kinetic theory of gases, the flow of a dilute gas is described by the Boltzmann 
equation (see e.g. [7, 9, 26]) 



where f(t,x,g) is the distribution function, and (t,x,£) 6 M. + x M. D x M D . Q(f,f) is the 
collision term with a quadratic expression given by 



(2.1) 




(2.2) 



where 



/*=/(o, /'=/(0: /: = /(&) 



(2.3) 
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and £' and are velocities after collision of two particles with original velocities £ and 
and with unit vector n 6 S^ 1 joining the centers of them, a is a constant equivalent 

to Na 2 , which keeps invariant when taking Boltzmann-Grad limit N — > oo, a — > (cf. 

[7, 9]). Here N is the number of particles while a is the diameter of each particle. 

However, such a collision term turns out to be too complicated for numerical simulation, 

so a variety of variants are raised to get it simplified. To some extent, the BGK operator 

[4] is the simplest one. It substitutes the collision term Q(f,f) by 

Qbgk(/) = -v(f - /m), (2.4) 
where v is the collision frequency, and /m is the local Maxwellian defined as 
„ , . . p(t,x) ( If — u(t, x)\ 2 \ 



[27r6{t,x)] D / 2 ^\ 20(t,x) 
It is related with / by 

p(t,x) = [ f(t,x,£)dt, 

p(t,x)u(t,x)= [ £f(t,x,£)dt (2.6) 

p(t,x)\u(t,x)\ 2 + Dp(t,x)9(t,x) = [ \t\ 2 f(t,x,Z)d£, 

where p, u and 6 can be viewed as macroscopic variables density, velocity and temperature, 
respectively. This model is much more easy to use in numerical methods. However, it 
suffers the disadvantage of being unable to predict the correct Prandtl number, which will 
be seen in the numerical tests. 

3 A numerical formation equivalent to Grad's moment method 

3.1 Discretization of the distribution function 

In order to solve the kinetic equations numerically, we first expands the distribution 
function into Hermite functions as in [12]: 

/(£) = E /«?M*)> ( 3J ) 



where a = (a±, ■ ■ ■ , old) is a D-dimensional multi- index, and 

V = —. 

Ve 

The basis functions Tlo,a are chosen as 

D 



(3.2) 



1 + 1 f \ 

H 6 , a (v) = J] -j=e-^He ad {v d )^ V (-^J , (3.3) 



where He ad is the Hermite polynomial defined by 



d" 



He n {x) = (-l) n exp(x 2 /2)^ r exp(- 2 ; 2 /2). (3.4) 
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Note that this is formally inconsistent with Grad's original expression, but will be more 
convenient for our deduction below 1 . The properties of Hermite polynomials can be found 
in many handbooks such as [1]. Some useful ones are listed below: 

1. Orthogonality: / He m (x)He n (x) exp(— x 2 /2) dx = m\V2Tv8 mjn ; 

Jm 

2. Recursion relation: He n+ \(x) = xHe n {x) — nHe n _i(x); 

3. Differential relation: He' n (x) = nHe n _\{x). 

It can be derived from the recursion relation and the differential relation that 

[He n (x) exp(- 2 ; 2 /2)] / = -He n+1 (x) exp(-x 2 /2). (3.5) 
The expansion (3.1) together with (2.6) yields 

D 

fo = P, f ei =0, ^/ 2ed = 0, i = l,---,D, (3.6) 

d=l 

where stands for the multi-index with the ith component 1 and all other components 0. 

It is known that (3.1) will result in an "infinite moment system". In order to make it 
numerically solvable, we choose a positive integer M ^ 2 and approximate (3.1) by 

/(£) ~ E ( 3 - 7 ) 

\a\^M 

Using Fm(u, 6) to denote the linear space spanned by all %0 )Q (v)'s, where \ a\ ^ M and v is 
defined by (3.2), then Fm{u, 0) is a finite dimensional subspace of L 2 (W N , exp(|t;| 2 /2) dv). 

Remark 1. Based on such an expansion, the Maxwellian in the BGK collision operator 
(2.4) can simply be expressed by 

/ M (£) = pH e ,o(v). (3.8) 

With (3.6), we have 

QBGK(/) = -f fcfreA v )- ( 3 - 9 ) 

K|a|^Af 

Remark 2. The definition of v (eq. (3.2)) implies that we take the mean velocity u as the 
"origin" and 9 as the scaling factor when discretizing the distribution function. Thus (3.6) 
holds. If (3.6) is violated, and we suppose u' £ M> D , 0' £ M + , and a distribution function 
/ is approximated as 

fig) « f«He<, a {v\ V = Lljf (3.10) 



1 Grad uses symbols as 3^j"i 2 ...j„ to denote basis functions. This symbol is equivalent to CH a , ce = 
ejj + e,: 2 + ■ • ■ + £;„ which is used here, where C is a constant factor. Inversely, a can be expressed as 

^ub ] with 

Qi ones, ct2 twos, • • ■ , and a_o D's in the subscript s?t6. Thus our expansion is actually the 
same as what Grad has done. 
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with / e . 's or fie nonzero, then the associated p, u and 9 can be calculated by substi- 
tuting (3.10) into (2.6). Since 

He (x) = l, Hei(x) = x, He 2 {x) = x 2 - 1, (3.11) 

employing the orthogonality of Hermite polynomials, the integrals on the right hand sides 
of (2.6) can be directly worked out as 

P = fo, 

P U = pu' + (fe d )d=l ,-,D, 

D (3-12) 

p\u\ 2 + Dpd = 2pu ■ v! - p\uf + Ytfh + 2 /2eJ- 

d=l 

Tolke, Krafczyk, Schulz and Rank [29] use it' = in their discretization. Additionally, the 
heat flux can be calculated by 

D 



\ I \i~ u\ 2 (^ - Uj) d£ = 30% + 2/ 3e . + V f ej+2ed , Vj = 1, • • • ,D. (3.13) 
2 J * D d=i 

3.2 Outline of the fractional step method 

In this subsection, our numerical scheme will be outlined. Suppose x G M. N and N ^ D. 
Th is a uniform rectangular mesh in M. N , with all grid lines parallel with the axes, and 
each cell is identified by an ^-dimensional multi-index /3. That is, for a fixed xq £ 
and Axj > 0, j = 1, • • • , N, 

T h = {Tp = x + \faAx!, {Pi + l)Axi] x • • • x [p N Ax N , (f3 N + l)Ax N ] : /3 e Z N }. (3.14) 

Using to approximate the average distribution function over the cell Tp at time t n , 

the Boltzmann equation (2.1) can be solved by a standard fractional step method: 

N At n 

1. Convection step: f^) = - £ ^[^(fl - ^_ i<y (01- 

3=1 3 

2. Production step: = f% +u (£) + Ai"Q h (/™ +1 *). 

In the convection step, the finite volume method is employed and F R , i_ is the numer- 
ical flux between cell and Tp +e . In the production step, Qh is a transform over 
Fm{u^ +1 *, 0p +1 *), and is considered as an approximation to Q(-, •). Here ■u^ +1 * and 0^ +1 * 
are the mean velocity and temperature corresponding to the distribution function 

The numerical flux can be chosen from the standard ones in the finite volume method, 
but in our framework, only central schemes are available since the characteristic factor- 
ization of the flux function is not available yet. A series of central schemes can be found 
in textbook [30], and we choose the HLL scheme [16] in our numerical experiments, which 
reads 



\R — \L , A^. < < Xj , 



(3.15) 
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Xj and are the fastest signal velocities arising from the solution of the Riemann prob- 
lem, which will be discussed in Section 3.3.2. For all /3 G Z N , given 

m) = E flavin = ce - ^)/(^) 1/2 , (3.16) 

where a = He n ,a, and u^,6^ are the mean velocity and temperature in cell /3, then 
£jfp (£) can be calculated according to the recursion relation of Hermite polynomials: 



' ' (3.17) 

= E ^ [w^^) + + «^,«,(«?) • 

Since |a + ej\ = M + 1 when |a| = M, £jf%(£) no longer exists in the space Fm(u^, 9p). 
Thus, we need an additional "projection step" to drag (3.17) back into Fm(u^,9^). This 
can be done by simply dropping the terms with \a + e,-| = M + 1, since when \a\ > M, 
T-L a {v) is orthogonal to Fm(u,0) with respect to the inner product 



(/,</) 



J f(v)g(v)ewP?p\dv. (3.18) 



However, the convection step is still uncompleted since it is nontrivial to add up two 
functions lying in Fm(u^,6^) and Fm ( u1 j3+ e ■ > @p+e ) respectively. This is a part of our 
major work and will be discussed in 3.3. 

As to the production step, the main job is to construct the numerical collision oper- 
ator Q h . This is also implemented by projecting Q{f^ +l \ f^ +1 *) into F M {u™ +u , fl™ +1 *). 
Precisely, Qh is defined as 

Qh(f)= E Q«Ho, a (v)> Vf£F M (u,0), (3.19) 

where 

Qa = C e , a f Q(f, f)(t)He, a (v) exp(M 2 /2) dv, (3.20) 

^e, a = : ; — • (3.21) 

ail--- a d l 

Further calculation requires the concrete forms of /' and fl in (2.3). For BGK model 
(2.4), the numerical collision operator has a simple explicit form (3.9). In this situation, 
all fa's can be decoupled, so the production step can be performed by solving each f a 
analytically. The scheme reads 

2. Production step (only for BGK model): 

Remark 3. For the time integration, we use a single step Euler scheme for both convection 
and production step. Actually, such formation can be smoothly generalized to Runge- 
Kutta and Strang splitting schemes. 
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Remark 4. For other collision terms, such as ES-BGK model [17] or Maxwell molecules [11], 
the expression of Qh can be much more complicated. For the ES-BGK model, it is always 
possible to get Q a 's by direct integration. For Maxwell molecules and linearized Boltzmann 
collision operator, the method in [34] can be employed to generate the numerical collision 
operator Qh- For simplicity, these models are not considered in this paper. 



3.3 Completion of the convection step 

Two points remain unclear for the convection step. One is that we need to find a 
way to add up two functions in different spaces Fm(ui,9x) and Fm(u2,02), so that it is 
applicable to calculate the numerical fluxes and to accumulate them to the solution at the 
last time step. And the other is the estimation of the characteristic velocities Xj and . 



3.3.1 Projection between two different spaces 

Assume f\ E Fm(u±,9i) and J% E Fm(u2,02)- Obviously, when u\ ^ 112 or 9\ 7^ 62, 
direct calculation of /1+/2 is inapplicable. Therefore, we want to find f\ E Fm(v>2, #2) such 
that fi is some approximation of /1 in Fm(u2,02)- In order to realize such transformation, 
we propose a fast projection method which has a time complexity of 0(M D ) below. 

First, let us consider the case of M = 00, and / E (iti , #1 ) H ^00(1^2, ^2)- Then, / 
has the following two representations 

/(£) = 2 h^He^ivi), t;i = (£- ui)/v^, (3.22) 

and 

/(£) = h«Ke 2 , a (v 2 ), v 2 = {£- U2)/VF 2 . (3.23) 



Suppose all /i, Q 's are known, and we want to solve all f2,a' s - Let 9 = y#i/#2 and 
w = (u\ — VL2)I\J~02- It is obvious that 

v 2 = 9vi + w, Hg 2 , a = 9 H+D n 9l)a . (3.24) 
Joining (3.22) and (3.23), we have 

Yl /i,«ft*i,a(«i) = Yl f2j lal+D U 9l A9v 1 +w). (3.25) 

Now we introduce an auxiliary function F(v,t), defined as 

F(v,t)= Y F a (r)[0-l)T + l]^ +D H ei ,a([0-l)r + l]v + Twy (3.26) 

aeN D 

which satisfies 

F a {0) = fi, a , Va E N D , and F(v,0) = F(v,l). (3.27) 

Comparing (3.26) (3.27) with (3.25), it can be found that for any a E N D , F a (l) is just 
f2,a which is to be solved. Moreover, if we suppose 

dF 

0^ = 0, VrE [0,1], (3.28) 
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then an infinite ordinary differential system of {F a {r)} Q& ^D can be obtained. 

dF 

The detailed calculation of — — can be found in Appendix B, and we only show the 

or 

final result here: 



a£N D I d=l 



e d + W d y%_F a - 



erf 



where 



R(r) 



1 



., S(t) = 1 - tR(t) 



(0-1)t + 1 (0-1)t + 1 

and the parameter of Hg lja is \{0 — l)r + l]v + rw. As required in (3.28), 



(3.29) 
(3.30) 



dr 



Q\RF a -2e d + WdV 6lF c 



a-e d 



d=l 



VaGN D , Vt€[0,1] (3.31) 



must hold. (3.31) is an infinite system, but for any M ^ 2, if we consider only a subsystem 
containing all equations with \a\ M, it is still closed. Therefore, in order to project a 
function / G Fco(**i, #i) H -Foo(w2, #2) into Fm(u 2 , #2)) it is only needed to solve (3.31) for 
all |q| ^ M. 

For an arbitrary function /(£) which is defined on MP , its projection to Fm{u,9) is 
defined as 



iw= y; c M / /(o^,a(^)ex P (i^i 2 /2)d^ 



1; = (£ - ii)/vft (3.32) 



See (3.21) for the definition of Ce iQ! . The following proposition provides an algorithm to 
project a function in Fm(u\,9\) to Fm(u2,&2)- 



Proposition 1. Suppose f G -Fm( u i>0i) can ^ e represented by 
/(£) = X] h,aHe u a{vi), vi = (£ - 

For some it2 G R #2 > 0, {-F q ,(t)}i q ,i^.m satisfies 

D 



(3.33) 



-^F Q = ^ S 2 ^iFF Q _ 2ed + w d V^X-eJ , Vr G [0, 1] 



d=l 



(3.34) 



1^(0) = fl, ai 

where S and R are given by (3.30), and w = (tti — i^Vv^- Let 

|a|<M 

T/ien G Fm(u2,02) and g(£) satisfies 



(3.35) 



Ml 



(3.36) 
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Proof. Let A = [—1, 1} D and 

PeuM = U^aM exp(|^i| 2 /2), Va G N D . (3.37) 
Thus {-P6»i,o(' l 'i)}aeN D forms a complete basis of L 2 (A). It follows that for an arbitrary 



set of {ba}\ a \<Mi there exists a function / G L 2 {A), such that 



/ /(«i)Ptf 1> a(«i) d«i = 6 Q , Va G N D , |a| < M. (3.38) 

When b a = f{^)Pe 1 a( v i) d«i, a zero extension of / onto R D gives that 
Jr d 

J /(ui)^i,a(«i) exp dwi = J f{^)'He 1 ,a( v i) ex P dvi ( 3 - 39 ) 

holds for all a G N D , |a| ^ M. Let <?(£) = /(i?i). Since g has a compact support on 

M D , we have g G .Foo(mi,0i) H F oo {u2,02)- Now (3.39) and the orthogonality of Hermite 
polynomials implies that if 

0(0 = J2 hc^OiAvi), (3.40) 

then g 1<a = f lj(X for all a G N D , \a\ < M. 
The preceding analysis shows that if 

9(Z) = Yl 92,an e2!a (v 2 ), (3.41) 

then §2, a = F a (l) f° r an a ^ f^^\ M ^ Employing the orthogonality of Hermite 
polynomials again, we can deduce that for any p(£) G Pm(K d ), 



/ 



p(£)/(£)d£ = / p(£M£)d£ = / P(£)5(£)d£. (3.42) 
Thus (3.36) is established. □ 

Based on this proposition, the projection requires to solve an ordinary differential 
system (3.34). Let row vectors -F , m (r) and F(t) be 

F m (r) = (F Q (r)) H=m , < m < M, (3.43) 

and 

F(t) = (F (r), Fi(t), ■ ■ ■ , Fm(t)). (3.44) 
Thus the ordinary differential equations (3.34) can be simplified as 

±F(t) = F(r)A(r), rG [0,1]. (3.45) 

Eq. (3.34) reveals that A(r) is an upper triangular matrix with vanished diagonal entries. 
Therefore, (3.45) can actually be solved by recursive integration. However, the direct 
integration will leads to 0(M 2D ) calculations, so we solve (3.45) by applying O(l) steps 
of Runge-Kutta numerical integration, which is unconditionally stable due to the special 
form of A(t). Since A(r) is sparse, each Runge-Kutta step costs only 0(M D ) calculations. 
Thus the whole projection has a time complexity of 0(M D ). 

Now let us return to the convection step. Using H^ to denote the operator that projects 
any function to the space FmM, 0q), then the convection step is described as following: 
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1. Convection step: 

(a) Apply the convection within Fm(u^,9^): 

N At" 



= m) - E^[( n ^; + |,)(C) - (n^;_ |ej )(0]; (3.46) 

(b) Use (3.12) to calculate the mean velocity it^** and temperature #^** for the 
distribution function /?**; 

(c) Let u n + u = uf*, = and = n™ +1 */^**- 

When implementing Step la, ITg is actually applied on each term of Fp±± e ., and the result 
fp** no longer satisfies (3.6). So the mean velocity and temperature need to be recalculated 
in Step lb. In Step lc, /£** is adjusted to such that (3.6) holds for Later on, 

in the production step, the mean velocity and temperature are not changed, so (3.6) still 
holds for fa + . The conservation of the convection step follows from proposition 1 and 
the conservative form of the finite volume scheme. To be specific, for any p(£) G Pm( 
we have 



E / p(o/r(o<* 



IV 



3=1 



At n 
Ax * 



TV 

E ^w-EE 



As, 



(3.47) 



[Using (3.36)] 



Thus quantities such as mass, total momentum and total energy are conservative. 

3.3.2 Estimation of the characteristic velocities 

In order to estimate A^ and A^ that are used in (3.15), we need to investigate into 
the expression of numerical flux njgi^ ± i e . carefully. Precisely, we should make sure the 
Riemann problem that such a numerical flux solves. In order to simplify the notation, we 
consider only the following form: 



Afn^/x) - X'fUj.jMj^h) + xfxf (U f2jl f 2 - h 

Xf-Xf 



(3.48) 



where f\ G Fm{u\,6i) and /a G ^Af (1*2^2)1 an d both satisfy (3.6). For / G Fm(u,9) 
which satisfies (3.6), H/ is the projection operator from F^u^O) to Fm(u,9), which 
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simply discards the terms with orders higher than M. II/ 2 a is the projection operator 
from Fm(u 2 ,9 2 ) to Fm(ui,6i). Then if we take 

h = fl f2 = f$ +ej , (3.49) 

F\ is exactly the same as H^F j3+ i e . in the case of < < Af. Similarly, let 



AflI/,/JI/:^/i) - A|n /2 (£,-/ 2 ) + Aj^Af (/ 2 - II/,/,/-,) 



Af-A- 



(3.50) 



Then F 2 is just n^i^_i if f x = ffi and f 2 = f$ in the case of Xf < < Aj*. Due to 
the similar forms of F\ and F 2 , only i*\ is considered below. 

The nature of Fi can be depicted with the help of the following proposition: 

Proposition 2. II/ 2 j x is invertible. 

Proof. Denote the projection operator from Fm{u\,9\) to Fm(u 2 ,9 2 ) as ILfaj^. We are 
going to prove that ILf 1 j 2 Hf 2i f 1 is the identity operator. Proposition 1 shows that for any 
feF M (u 2 ,e 2 ) andpG P A/ (i D ), 

I K0/(0^= I P(0(n f2l/l /)(0d^= I P ^)(u fl j 2 u f2jl m)dc (3.5i) 

jrd Jm. d Jr d 

That is, 

p(0[(/ - n /lj2 n /2i/l )/]^) d^ = o. (3.52) 



/ 



Choosing = T-L0 l>a (vi) exp(|i>i| 2 /2) for a G N D , |a| ^ M respectively, and making 
use of the orthogonality of Hermite polynomials, it follows that 

(J-n /li/2 n /2)/l )/ = 0, Vf £ F M (u 2 ,0 2 ). (3.53) 

Similarly, it can be proved that n^j^II/^/,, is also the identity operator. Thus II/ 2) / 1 is 
invertible. □ 

Now let us turn back to the numerical flux (3.48). Let f 2 = Hf 2 j 1 f 2 - Based on 
proposition 2, we rewrite (3.48) as 

F Afii / . / ii / ^ f ii /: :,/ l i - xfn hJl u h ^uj 2 \ fi f 2 ) + \f\*{} 2 - a) 
1 = x R - \ L ' ^ ^ 

where n^ 1 j 1 is the projection from Fm(ui,9i) to itself, which is actually the identity 
operator. Now it is clear that the corresponding Riemann problem of F\ is 

d f d 

i + ^[n/,/ 1 n / (e i n / -} i )/] = o ) 

|/(0,x) = / 1) *<0, ^ 3 - 55 ) 
\/(0,x) = / 2) x>0. 

Here / always lies in Fm(u\,9i), and the meanings of II f ^ and II j have been changed a 
little. Suppose u and 9 are the mean velocity and temperature associated with /, whose 
explicit expressions can be obtained from (3.12). Then 11/ a is defined as the projection 
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operator from Fm(u,6) to Fm(ui,6i), and Hf is denned as the projection operator form 
F M+ i(u,e) to F M {u,6). 

The characteristic velocities of Riemann problem (3.55) seem to be difficult to obtain. 
Therefore, in order to give an estimation of Xj and A^, we choose a fixed distribution 
function /* E Fm{ui,9\) that lies "between" f\ and J2, and linearize (3.55) as 

% + n / , /l n / .(f J n- fi |)=o. (3.56) 

Thus, we only need to estimate the eigenvalues of II/* j x TLf*£jHJ* - , which is an operator 
on Fm{u\,6i). Since II/*,/! is linear and invertible, the problem can be further simplified 
as the estimation of eigenvalues of Hf*£j, which is an operator on Fm(u*,8*). Taking 
v* = (£ — u*)/\fd*, we have 

n f . =U f * (u* + v* V¥) = u*I + V¥U f . v* . (3.57) 

For the eigenvalues of Uf*Vj, we have the following proposition: 

Proposition 3. The eigenvalues ofHf*v* are formed by all the zeros of He m +i(x), m = 
0, -- ,M. 

Proof. Suppose m E {0, • • • , M}, and all zeros of He m+ i{x) are denoted as xq,--- ,x m . 
Note that all x,'s are real and different (see e.g. [25]), so we can assume that 

x < ■■■ < x m . (3.58) 

For any i £ {0, • • • , m}, there exists a unique polynomial Pi, m (x) G P m (x) that satisfies 
Pi,m{xk) = <5jfc, k = 0, • • • , m. Let a € N D satisfy |q| = M — m and ctj = 0. We are going 
to prove 

^■f*[vjPi,m(Vj)'He*,a(v*)] = x i pi > m{v*j)'He*, OL {v*). (3.59) 
The definition of n f* shows 

Uf*[VjPi tm (Vj)H9*, a {v*)] = V*Pi tm (Vj)H9* t a{v*) - Cj, a He m+1 {Vj)n9*, a (v*), (3.60) 

where Cj i0l is a properly selected constant such that (3.60) lies in Fm(u* ,9*). Thus, for 
any k E {0, • • • , m} and v* k satisfying v^. ■ = xy., we have 

Ilf*[v*Pi,m(v*)ne*, a (v*)}\ vlt=v ^ = x k 5ikHe*,a{v* k ). (3.61) 

Then (3.59) holds due to the uniqueness of Pi^n- 

It remains to prove that Ilj.^j has no other eigenvalues. Let us count how many 
eigenvectors are included in the form of Pi ! m(Vj)'He* :a (v*). Consider an arbitrary a E N D 
and \a\ ^ M. Let 

i = ctj, a = a — iej, m = M — \a\. (3.62) 

Obviously i, a and m satisfy ^ i ^ M and ctj =0. Thus each a uniquely determines 
an eigenvector. However, the number of such d's are equal to the dimension of space 
Fm(u* ,0*), so Iif*Vj has no other eigenvectors, thus no other eigenvalues, either. □ 
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According to Proposition 3, the smallest and largest eigenvalues of Tlf*v* are the 
smallest and largest zeros of HeM+i{%), denoted by xq and xm- (3.57) shows that the 
smallest and largest eigenvalues of LLf*£T are u*j + xoV&* and u* + xmvO* respectively. 
Since /* lies "between" f\ and f 2 , we use 

Xj = minjitij + xoy/(h, u 2 j + xq\/62}, (3.63) 
A^ = maxjnij + x M \fOi, u 2 j + im\/^} (3.64) 

while computing numerical fluxes. In our implementation, we use the subroutine in AL- 
GLIB [6] to calculate the roots of HeM+i, and Xj and X^ are also used in the CFL 
condition to determine the time step length. 

Remark 5. The numerical method described in this section is only of the first order. In 
order to extend it to higher order schemes, reconstruction techniques need to be added 
to the finite volume scheme. Since addition and subtraction between two distribution 
functions are already available, it is only needed to determine a proper "slope" , which can 
probably be done with the help of the standard slope limters used in the normal finite 
volume schemes. 

3.4 Relation with the Grad-type moment method and the LBE model 

The relation and difference between the Grad-type moment method [12, 14] and the 
LBE (lattice Boltzmann equations) model [28] are summarized in [24], where both mod- 
els are considered to be some approximation to the Boltzmann equation by an Hermite 
polynomial expansion. The expansion (3.1) and truncation (3.7) are exactly the same as 
what have been done by Grad [12], which means the method described above is actually 
solving the Grad-type moment equations. For M = 3, 4, 5 and D = 3, it corresponds 
to the 20, 35, 56-moment equations which take the complete Mth order moments. How- 
ever, systems such as 13, 26, 45-moment equations are not included. Those complete Mth 
order moment equations are popular in extended thermodynamics, see e.g. [22, 33, 3]. 
A software called ETxx [2] is developed by J. Au, H. Struchtrup and M. Torrilhon to 
generate equations for arbitrary order of moments, however, explicit moment equations, 
when written in conservative form, require 0(M 2D ) calculations for the flux function, and 
this is reduced to 0(M D ) in our numerical formation. 

Now that our numerical strategy is equivalent to Grad's moment method, one can refer 
to [24] for the precise difference between our method and the LBE model. According to 
[24], the linearized equation (3.56) can be interpreted as the LBE model, which places the 
center of the lattice at u* . Thus such linearization is reasonable. 

4 Regularization of the moment method 

The main drawback of Grad's moment method is that its hyperbolicity yields unphysi- 
cal subshocks [13]. The behavior in the case of high order moment equations can be found 
in [33, 3]. [14] provided a way to regularize Grad's moment equations, and it was further 
studied in [27, 36] as R13 equations. Now we follow the regularization technique in [27] 
and regularize the numerical method introduced in Section 3 in exactly the same way. 
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4.1 Chapman-Enskog expansion around the truncated distribution 

Independent of Grad's method, Chapman-Enskog expansion [8, 10] is another im- 
portant method for deriving equations of macroscopic variables. Following the generic 
procedure of Chapman-Enskog expansion, a scaling parameter e is introduced on the right 
hand side of the Boltzmann equation. However, according to [27, 26], only the high order 
part instead of the whole collision term is scaled: 



df 



1 



^ + ^-VJ = Q(/°, /°) + -[Q(f, f) - Q(f°, f% 
where f° is a truncation of the distribution function (3.1) defined by 

\a\^M 



(4.1) 



(4.2) 



The scaled part Q(f, f)—Q(f°, f°) will be denoted as Q below. Let us apply the Chapman- 
Enskog expansion around / , i.e. expand / by 



f = f° + ef 1 + e 2 f 2 + --- , 



(4.3) 



and we require that the truncation at any term of this expansion keeps the same values 
of all the moments with orders less than or equal to M. Suppose the scaled part of the 
collision term Q has a corresponding expansion 



Q = eQ 1 + e 2 Q 2 + 



(4.4) 



It is reasonable to assume that Q has no zeroth order term since it has been taken away 
from Q(f, /). Match the zeroth order term on both sides of (4.1), and we have 



Q 1 



dt 



+ ^V x f°-Q(f°J ). 



(4.5) 



For any multi-index a with \a\ > M, multiplying ~Hq a ( v ) ex P(M 2 /2) on both sides and 
then integrating the whole equality over R D with respect to v, since 



9/0 exp(M 2 /2) = d[/° ex P(H 2 / 2 )] _ /0 <9[exp(M 2 /2)] 



dt 



dt 

d[/°exp(H 2 /2)] 
dt 



+ 



u 



dt 





Ve dt 



u 



(4.6) 



/°exp(M 2 /2), 



the orthogonality of Hermite polynomials leads to 
^ ia (^)g 1 exp(|^| 2 /2)d^ 



^,a(^)G(^/ )exp(| U | 2 /2)d^ 



^,a^)Q(/ ,/ )exp(|«|72)d«, 



where 



J_d_ fv_\ 



an 



(4.7) 



(4.8) 



Now the concrete form of the production term is required for further calculation. Still, 
we adopt the simplest BGK model, which gives 



Qbgk(/°) 



-K/°-/m), 



Q 1 = -vf l - 



(4.9) 
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Thus the second term on the right hand side of (4.7) also vanishes, and (4.7) becomes 

/ ^,aH/ 1 exp(|^| 2 /2)d^ = -i / n e , a (v)G(a°)ew(N 2 /2)dv. (4.10) 



This shows that f 1 can be represented as 



V 

\a\>M 

where the coefficient [G(£f°)] a is the corresponding coefficient of G(£/ )'s expansion in 

Foo{u,e). 

At last, we set e = 1 and approximate / by / « /° + / . Since f 1 can be obtained 
from /° which is of finite dimension, when substituting such / into the Boltzmann-BGK 
equation, a closed system can be obtained without more truncations. For D = 3 and 
M = 3,4, 5, the R20, R35, R56 equations for the BGK model can be obtained. 

4.2 The numerical method 

In this subsection, we restrict our focus on the BGK model. As in Section 3, only 
f° is stored at each time step. Note that for the BGK model, f° and f 1 are decoupled 
in the production step, which can be implemented exactly the same as that in Section 3. 
Therefore, we concentrate only on the convection step below. 

Consider the original form of the numerical flux (3.15), and now fS is recognized as 

f/3 +fp> where 

ffEF M {u%ej), n^ = 0. (4.12) 

The latter equation follows from (4.11), (3.32) and the orthogonality of the Hermite poly- 
nomials. Similarly, we have 

n ^± ej =°> V/3 £ II 1 . (4.13) 

Since only is desired after the convection step, the form of (3.46) is still adoptable. 

Therefore, we can mimic (3.48) and write the numerical flux as 

„ ffljftfo/jj) ~ Ajdlr,/,!!^^,/^; + Aj Afill;. ^ - fj) 

Xf-Xf 

+ Afl'/ i^/i'l " A' lI / ,. /l llr ! ii ; ./^) (4.14) 



Xf-Xf 



Fn + F12, 



where /i = + /f\ ji = + f\ , and Tl^fl = n^/j = 0. As described in Section 3, Fn 
is just the flux for Grad's moment equations. As to F±2, (3.17) and (4.11) show 

= — ( ai J + 1 )[ Gr (€/°)] tt +e,-^ i ,a(«i), * = 1,2. (4.15) 

\a\=M 

It is easy to see that in the expansions of f\ and /2, only the coefficients with \a\ = M + 1 
have effect on the numerical flux F±. We can also find that when M ^ 3, F12 in (4.14) has 
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actually no contribution to the velocity u and temperature 9 for the next time step, since 
(4.15) reveals that the Grad's expansion of F\2 contains only terms with orders higher 
than or equal to M, but (3.12) tells that u and are only relevant with the coefficients 
with orders less than or equal to 2. That is to say, for all /3's, u^ +1 and can be 

solved using the method introduced in Section 3, without adding the "regularizing part of 
numerical flux" F\2- Thus, the time derivative in (4.8) can be explicitly approximated by 




(4.16) 



Now, in order to approximate [G(£f°)] a , it is only needed to approximate 



[V* • (£/%, ol e N D , 
For a fixed point Xq £ R^, we have 

d 



a 



M + 1. 



(4.17) 



TV 

[v, • (af)Ux Q ) = 

N 

=£ 

N 

= £ 

N 
3=1 



&7°) 



(^o) 



d dxj 



d 
dxj 



&f o )n 0o , a dv 



(v 



x=xq 



(4.18) 



^0, « 



x=xo 



X = XQ 



where Ce 0>Q is defined by (3.21), uq and 9q are the mean velocity and temperature at 
point £Co, and n uOi o is the projection operator to the space Fm+i{uq,0q). Now consider 
the discrete circumstance. For each /3 6 Z N , ^u^,d n (CjfJj± e ) can De obtained according 
to (3.17) and (3.34). Without confusion, the superscript "n" will be omitted below, and 
n U fl,0o is simplified as IL3. Mimicing the method in [32], suppose 



di 



Ax, 



Ace; 



J / o 



(4.19) 
(4.20) 



and then we can reconstruct the spatial partial derivative by a central difference in the 
smooth case: 



di + d.2 



P 2 

or the van Leer reconstruction in the discontinuous case: 



s Xj n^r) 



van Leer ^ l^cfo + l^l^l 
/3 |<il| + |d2| 



(4.21) 



(4.22) 
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Thus we have 



i ( N 

|a|=M+l I j=l 




J_,.n+\ n \ 



where (3.17) has been incorporated into the expression of G(£f°). 
From (4.15), we can conclude that 



(4.23) 



n /3 (^/ /3 1 ) = n /3 fen^ 1 ), 



j 



i, 



(4.24) 



which appears twice on the right hand side of (4.14). This expression can be computed 
cheaply since ILpfjj has only 0(M D ^ 1 ) non-zero coefficients. The last term in (4.14) 
requires us to compute 

n^n^e.^-/^.), (4.25) 

i.e. to apply another projection to the result of (4.24). Note that F a does not appear on 
the right hand side in (3.34), and for all a 6 N D , \a\ ^ M, the coefficients of Tie^^ m the 
expansion of (4.24) are zero, so all coefficients do not change after projection, saying if 



^I3±ej{tj f}± ej ) ~ ^2 C ±,<*H-8p± a(V/3±e 



(4.26) 



then 



\a\=M 

n^n /3±e .(^/ / ^ ±e .)= J2 C±,aHe p , a (vp)- (4-27) 

\a\=M 

e d +i 

This is the reason why a coefficient J 2 is multiplied in the definition of basis functions 
(3.3). Until now, the calculation of numerical fluxes for the regularized moment equations 
is thoroughly clarified. 

Remark 6. In (4.19) and (4.20), ^/3{Cjfjs± ej ) needs to be computed, while in the convection 

term of Grad's moment equations, or rather, F\\ in (4.14), Rp(fijfp± e .) also needs to be 

computed. Since Hp = TLpILp, and the second Hp is trivial, these two projections can be 
combined into one. 

Remark 7. As is known, the regularized moment equations contain second order derivative 
terms, so the CFL condition for the method above is 



N \ ■ 

A J,max 

^ Axi 

7=1 3 



1 + 



Ax a 



< 1, 



(4.28) 



where Aj jmax is the maximum of |A^| and \X^\ on all cells, and ^ max is the maximal collision 
frequency. Note that for this regularized model, when calculating \j and \^ using (3.63) 
and (3.64), the roots of HeM+i(x) should be replaced by the roots of HeM+2(x) since we 
use the (M + l)-th order moments of f 1 . Eq. (4.28) leads to a relatively small time step 
length. The time step length can be enlarged following the methods in [39, 32, 35], which 
is not yet implemented in our program. 
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Remark 8. In our implementation, the classical fourth-order Runge-Kutta method is used 
to solve (3.34). Since the solution of regularized moment equations is generally smooth, 
most of the time, only one Runge-Kutta step is able to provide enough accuracy for such 
local projections. In the case of sharp initial values, some more steps are performed. 
However, such a situation only appears at the very beginning of the calculation. 

4.3 Outline of the algorithm 

As a summarization, our numerical method for the regularized moment equations is 
outlined as below: 

1. Let n = and set the initial state fa' for all /3's. 

2. Determine the time step length according to the CFL condition (4.28). 

3. Apply the convection step in Section 3. 

4. Obtain flg/^' 1 as in Section 4.2. 

5. Add the "regularizing part of numerical flux" to fg + '°. 

6. Apply the production step at the end of Section 3.2. 

7. Let n <— n + 1, and return to step 2. 

4.4 Boundary conditions 

Currently, boundary conditions are not available in our numerical scheme. The bound- 
ary condition is always a delicate issue in the deduction of the macroscopic equations. 
Generally speaking, the kinetic boundary condition introduced by Maxwell [19] is ex- 
pected to be added to the algorithm. As in [12] and [37], half-space integration needs to 
be performed during the construction of the distribution function in the ghost cells. This 
is possible to obtain for the discrete distributions /° + il u e/ 1 , thanks to the recursion 
relation of the Hermite polynomials, but such integration requires a subroutine with a 
time complexity of 0(M 2D ), which results in much more computational time on the cells 
next to the wall. The details are still in preparation. 

5 Numerical examples 

In this section, ID and 2D numerical examples of our method for the regularized 
moment equations are presented. In all these tests, the global Knudsen number is denoted 
as Kn, and the collision frequency v (see (2.4)) is substituted by p(t,x)/Kn. The CFL 
number is always 0.8. We use the POSIX multi-threaded technique in our simulation, and 
at most 8 CPU cores are used. 

5.1 One dimensional case 

Two ID examples are studied as follows. Since the boundary conditions are currently 
out of our consideration, only free or periodic boundary conditions are used in the following 
examples. In this section, some numerical solutions of the Boltzmann-BGK equation are 
provided, which are obtained according to the algorithm described in [20]. 

5.1.1 Shock tube test 

The shock tube test has been investigated in many works due to its fundamental role in 
characterizing the hyperbolicity of equations. For Grad's 13-moment system, it is shown 
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in [31] that unphysical subshocks can be found. However, in [32], the numerical result of 
the Riemann problem shows that R13 equations are able to capture these waves correctly. 
Here we first repeat this test in [32] to obtain similar results. 
As in [32], the initial conditions are 



p(0,x) 




p(0,x) 




(5.1) 



where the pressure p equals to p8. The velocity is zero everywhere, and the fluid is 
in equilibrium everywhere. The computational domain is [—1,1]. In order to make a 
comparison, we set Kn = 0.02 and compute this model until t = 0.3 as in [32]. The result 
of R20 equations (M = 3) is given in Figure 1. Compared with the result in [32], the plot 
of density agrees with that in [32] very well. The plot of heat flux has a similar shape 
with that in [32], but they differ in magnitude. This is due to the different models in the 
collision operator. The BGK model fails to predict the correct Prandtl number, which 
results in the incorrectness of heat flux. 




Figure 1: R20 results for the shock tube test with Kn = 0.02. The thick line with the left 
y-axis is the plot of density, while the thin line with the right y-axis is the plot of heat 
flux. 1000 grids are used for calculation. 



Since Kn is small, almost the same results are produced by M > 3 and no subshocks 
are found. This illustrates the capability of our method in capturing physical waves. Note 
that the peak and valley of the heat flux lies in the "discontinuities" of the density, which 
indicates the non-equilibrium. The negative part (x < 0) has smaller heat flux due to the 
larger density or collision term. 

Now we set Kn = 0.5 to investigate the numerical behavior of our method in the case 
of greater Knudsen number. The curves of density and temperature for M = 3 to 8, as 
well as the numerical solutions of the Boltzmann-BGK equation, are plotted in Figure 
2. With this Knudsen number, the hyperbolicity of the regularized moment equations 
clearly turns to be dominant. The solutions in Figure 2 exhibit a similar behavior with 
those in [33]. In the region of x < 0, all results are similar because of the high density 
behind the initial shock. In front of the initial shock, both the density and temperature 
are converging to the BGK solution, although the convergence rate is much slower. 
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(a) R20, M = 3 (b) R35, M = 4 




-1 -0.5 0.5 1 1.5' -1 -0.5 0.5 1 1,5 



(e) R120, M = 7 (f) R165, M = 8 

Figure 2: Results for the shock tube test with Kn = 0.5. The thick line with the left y-axis 
is the plot of density, while the thin line with the right y-axis is the plot of temperature. 
The dashed lines are the numerical solution of the BGK equation. 1250 grids are used for 
calculation. 
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5.1.2 A test with smooth initial values 



This example is again from [32]. The initial conditions are 

p(0, x) = 2 + - cos(-7nr), u(0, x) = ^1 + - sin(7ra;), - sin("7rx), 0^ , (5-2) 

and the fluid is in equilibrium everywhere with p(0,x) = 1. Periodic boundary condition 
is used and the computational domain is the interval [—1,1]. In order to validate our 
method, we use Kn = 0.01,0.1,0.5 in our numerical computation exactly as in [32]. The 
end time is t = 0.4. 

The results for different Knudsen numbers and different moment equations are plotted 
in Figure 3. All tests are computed using 1000 grids. The results in the first column are 
almost identical, which indicates the correct behavior of our method in the dense limit. 
The R20 equations, which should be the closest to the R13 system, produce similar results 
as those of R13 reported in [32]. The temperature plots in the first row can be used to 
make comparison. 

The new results are presented in the second and third columns, where the numerical 
solutions for high-order moment equations are listed. For Kn = 0.1, the R20 result shows 
an incorrect profile of density. With increasing M, both the density and the temperature 
tend to converge. For Kn = 0.5, the R20 equations provide completely wrong structures, 
although smooth initial values are used. Results for even larger moment systems are 
plotted in Figure 4. In this plot, the satisfying temperature plot is obtained when M = 7, 
but the density plots behave similarly as those in Figure 2. The curves for odd and 
even order of moments hold different profiles, and they are toddling close to each other 
gradually. With M been increased up to 11, the density curve eventually exhibits a 
satisfying convergence. The phenomenon illustrates the necessity of large moment systems 
in the microcase. 

5.2 Two dimensional case 

Two 2D examples are investigated in our numerical simulation. Both examples use 
uniform grids in the spatial discretization. Though much more computational cost are 
needed for 2D problems, the equations with up to 84 moments are considered. 

5.2.1 Shock-bubble interaction 

In this section, the shock-bubble problem tested in [32] is repeated. The initial state 
contains a shock wave at x = —1.0 travelling with Mach number Mq = 2.0 into an 
equilibrium area with (p,u,9) = (1,0, 1). A bubble is in front of the shock with density 
profile 

p(0,x) = 1 + 1.5exp(-16|a; - x \ 2 ), (5.3) 

where xq = (0.5, 0) T , and constant pressure p = 1. The shock wave has a fully developed 
structure instead of a discontinuity. Thus a pre-computation of the shock profile is needed. 
The initial density surfaces for Kn = 0.05 and Kn = 0.1 are shown in Figure 5. A uniform 
mesh with 1000 x 400 grids is used in our numerical simulation. 
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Figure 3: The numerical results for problem (5.2). The thick line with the left y-&xis is 
the plot of density while the thin line with the right y-axis is the plot of temperature. 



The shock structure can be obtained by solving a ID Riemann problem constructed 
according to the Rankine-Hugoniot condition. The left state is 

4M n 2 / /5 Ml + 3 „ \ T 5M n 2 - 1 



and the right state is 

p r = l, u r = ^-^M ,0,oj , p r = l. (5.5) 

Both states are in equilibrium. After a sufficiently long time, a stationary shock will form. 
It is quite convenient to transform a stationary shock to an unstable one in our numerical 
framework. Suppose a ID steady shock is presented by 

= E f<*WHei*)A v ( x ))> *>(*) = - ^=r , ^R. (5-6) 
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(a) R120, M — 7 (b) R165, M = 8 




(e) R364, M = 11 (f) R455, M = 12 

Figure 4: The numerical results for problem (5.2) in the case of Kn = 0.5. The thick line 
with the left y-axis is the plot of density while the thin line with the right y-axis is the 
plot of temperature. The dashed lines are the numerical solution of the BGK equation. 

and it satisfies (3.6). Then, for an arbitrary velocity s = (s,0, 0) T , let 

£ - u'(x) 

u'(x) = u(x) + s, and v'(x) = , . (5.7) 

y/9{x) 
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(a) Kn = 0.05 



(b) Kn = 0.1 



Figure 5: The initial plot of density 



Substituting u and v by u' and v' in (5.6) and keeping all the coefficients unchanged, then 
(5.6) becomes an unsteady shock travelling with speed s. Let s = (y/bMo/3, 0, 0) r , and 
then the desired shock wave can be generated. The initial values of 9 for Kn = 0.01, 0.05 
and 0.1 at y = are plotted in Figure 6. 



2 - 



1.5 - 



1 - 



0.5 - 




Kn.0.01 

Kn.0.05 

Kn.0.1 



-2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 3 



Figure 6: The initial values of temperature at y = 

This example is aimed at the validation of our algorithm in the 2D case. To make 
comparison with the results in [32], Kn = 0.01, Kn = 0.05 and Kn = 0.1 are considered 
and only the R20 equations are simulated since it is the closest moment system to R13. 
Results for the dense case Kn = 0.01 at t = 0.8 are shown in Figure 7. Comparing with the 
results in [32], the profile exhibits a qualitatively agreement while the peak value between 
x = 0.5 and x = 1 disagrees. This is believed to be caused by the highly dissipative 
numerical flux without gradient reconstruction in our implementation. For Kn = 0.05, 
our R20 results of density and temperature (Figure 8) are much closer to those presented 
in [32]. But again, the heat fluxes show the same profile with different magnitude, owing 
to the BGK model used here. The whole structure after interaction with the bubble is 
drawn in Figure 9, with a good agreement with the former results. 

For Kn = 0.1, we know from Figure 3 that R20 results deviate from the BGK solution 
slightly, so some deviation between R13 and R20 results is reasonable. Our R20 results 
are plotted in Figure 10. A comparison with R13 results in [32] shows that both R13 and 
R20 equations are able to give correct structures of density and temperature, while NSF 
is not. 
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(a) Density plot (b) Temperature plot 

Figure 7: R20 results of shock-bubble interaction for Kn = 0.01 at t = 0.8 and y = 




(a) Density (b) Temperature 

Figure 9: R20 results of shock-bubble interaction for Kn = 0.05 at t = 0.9 





(a) Density plot (b) Temperature plot 

Figure 10: R20 results of shock-bubble interaction for Kn = 0.1 at t = 0.8 and y = 



5.2.2 An example with three-dimensional velocity 

In all the numerical examples above, the z-component of velocity is always zero. Now 
we consider an example with three-dimensional velocity with initial conditions as 



p(0, x) = 2 + - cos(7ra;i) + - sin(7rx 2 ), p(0, x) 
/ 1 + \ sin(7rxi) + \ C0S(7TX2) 

u(0, x) = i sin(7rxi) + ^ cos(7rx2) 
\ \ sin(-7rxi) + i cos(-7ra;2) 



1. 



(5.8) 



The fluid is in equilibrium over the whole computational domain [—1,1] x [—1,1] together 
with periodic boundary condition. 

For this example, the simulations of the R20 and R84 equations with Kn = 0.1 are 
carried out. The numerical solution of the R84 equations on meshes with different sizes 
are compared to check the spatial convergence order of our scheme. In the case of no exact 
solution being available, we take the numerical result on a mesh with 500 x 500 grids as 
the reference solution. Other results are computed on meshes with N x x N x grids, where 
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N x = 10 up to 200. HLL flux without gradient reconstruction is used in our finite volume 
scheme, so the convergence rate is expected to be the first order. 

The numerical solution is shown in Figure 11 and Figure 12. For density and temper- 
ature, R20 and R84 results are almost identical at both t = 0.2 and t = 0.4. However, 
observable deviation appears in the vertical heat flux in both Figure 11 and Figure 12. The 
nontrivial vertical heat flux declares the capacity of our method to simulate 3D non- 
equilibrium processes. The L 1 errors of the solutions on different meshes are illustrated 
in Figure 13, where E(-) is calculated by 

jvi ref) A^ ref) 

EW = log 10 £ £ Axi rcf) Axr f) |^ {num) (^ ( f ) - V (rcf) K ( f )l- (5-9) 

i=l j=l 

Here all symbols with superscript "(ref)" stand for the corresponding quantities in the 
reference solution, i.e. the solution on the 500 x 500 mesh, and the symbol with superscript 
"(num)" is the solution on the coarse mesh, which is considered as piecewise constant. 
Obviously, first order convergence rate is achieved. 



Density p Temperature 6 Vertical heat flux q% 




Figure 11: Numerical solution of problem (5.8) at t = 0.2 and y = 0. 



6 Concluding remarks 

A uniform method to solve the regularized moment equations for arbitrary order is pro- 
posed. This is the first time that the method for arbitrary order regularized moment equa- 
tions is developed, and the moment method for large systems is applied to two-dimensional 
problems. We are now devoting our efforts to the mesh adaptation and parallelization of 
the algorithm to improve the computational efficiency so that the proposed method can 
be applied to practical applications. 
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Density p 



Temperature 9 



Vertical heat flux q% 





Figure 13: The error plots for problem (5.8). The x-axis is the logarithm of N x , and the 
y-axis is the logarithm of the norm of the L 1 error. 
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Appendix 

A Collection of the mathematical symbols 

Since lots of mathematical symbols are used in this paper, in order to provide conve- 
nience to the readers, we list some of them here. 
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Q(f,f) The Boltzmann collision operator 

Qbgk(/) The BGK collision operator 

/m The Maxwellian distribution 

He, a The basis functions for Grad's expansion 

He n {x) The Hermite polynomials 

Fm(u, 9) The finite dimensional space spanned by %0 )a ((£ — u)/y/6), \a\ ^ M 

fa The discrete distribution function on the cell indexed by (3 at time t n 
The mean velocity and temperature on the cell indexed by (3 at time t n 



,n on 



Fp + i e . The numerical flux between cells indexed by /3 and /3 + ej 

Xj 1 , X^ The fastest signal velocities travelling in the direction of — xj and xj 

Qh(') The discrete collision operator 

a Equivalent to r Hff^,a 

C e , a See (3.21) 

n u ,6»/ The function generated by projecting / into Fm(u,9) 

Abbreviat ion of Yl u n On 

Hf The projection operator from F^it, 6) to Fm(u,0), where u and 6 is 

the mean velocity and temperature of the distribution function / 

Ilftjj The projection operator from Fm{u\,6\) to Fm(u2, 2), where ui and 
Oi are the mean velocity and temperature of the distribution function 
fi, * = 1,2 

/° A truncation of the distribution (3.1), defined by (4.2) 

f k The fcth order term in the Chapman-Enskog expansion 

(■) a The coefficient indexed by a in the expansion of the parameter function 

H u ,o The projection operator to the space Fm+i(u,9) 

Up Abbreviation of Il^n^" 



B Calculation of the partial derivative dF/dr 

The calculation of the temporal partial derivative of (3.26) is performed here. Define 
A(v,w,t) and B(v,w) as 

OA 

A(v,w,t) = [{9- 1)t + 1]v + wt, B(v,w) = — = {9 - l)v + w. (B.l) 



It follows from (3.5) that 
d_ 



\A(v,w,t)] 2 

He m (A{v,w,T))exp 1 



B(v, w)He m+1 (A(v, w, r)) exp 



2 

[A(v,w,r)f 



(B.2) 

2 

With the definition of R(t) and S(t) (3.30), B(v,w) can be related with A(v,w,t) by 



B(v,w) = R(t)A(v,w,t) +wS(t). (B.3) 
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Substituting (B.3) into (B.2), and employing the recursion relation of Hermite polynomials, 
we have 

^[He m (A)eM-A 2 /2)] ^ 

= - [RHe m+2 (A) + wSHe m+l (A) + (m + \)RHe m {A)} exp(-,4 2 /2), 

where all parameters of A, S and R are omitted. With (B.4), the partial derivative of 
Hq 1jOC (A(v, w, t)) can be naturally obtained: 



— H &1 , a (A) =-J2 [e 1 R'He 1 , a +2e d { A ) + W d ^SH ei , a +e d ( A ) + ( a d + ^)R^,a{ A ) 

(B.5) 



d=i 



where 

A = A(u, w, t) = [(§ - l)r + l]v + tw = [A(ux, w±, r), ■■■ , A(u D , wd, t)] t . 
Since 



we finally get 



F(v,r) = F Q (r)5(r)-(l Q l +D )^ lia (A(^, W ,T)), 

a£N D I 
D 

d=l 

+ (\a\+D)(l-0)SHg lia 

aSN D I d=l 



(B.6) 
(B.7) 



(B. 



where the parameter of Tie 1 , a , i-e. A(u,w,t), is also omitted, and for a with negative 
components, F Q is taken to be zero. 
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